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Abstract 



'-^J \ In this work a three-states cellular automaton proposed to describe part of a biological immune 

system |y| is revisited. We obtain the dynamic critical exponent z of the model by means of a recent 
technique that mixes different initial conditions. Moreover, by using two distinct approaches, we have 
, also calculated the global persistence exponent 9 g , related to the probability that the order parameter 
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CN) ' of the model does not change its sign up to time t [P(t) oc t ° B ]. 
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I. INTRODUCTION 



Cellular automata are statistical-mechanical models which present a complex behavior, de- 
spite their relatively simple dynamic rules. In general these models are defined in ci-dimensional 
lattices with linear dimension L, in which each site % (1 < i < L d ) is occupied by a variable Sj. 
The dynamic rules of a typical cellular automaton have two basic characteristics: Firstly, the 
update of a variable Sj depends only on its neighborhood (short-range interactions) . Secondly, 
given a configuration of the system at instant t, at instant t + 1 all variables (1 < i < L d ) are 
updated at once. The use of cellular automata in order to mimic biological immune systems has 
improved the understandingabout the microscopic mechanisms that lead to the macroscopic 
behavior of these systems For example, cellular automata are as usefull to describe the 
T-helper cells response under parasitic infections Q, 0] as they are to model the course of the 
Human Immunodeficiency Virus (HIV) infection |^|. 

Brass et al. divised a simple cellular automaton model of Th cell interactions to mimic the 
immune system of mice exposed to parasitic infections . The model is defined in a cubic 
(d = 3) lattice in which each site is occupied by one of three different Th cell types: ThO, ThI 
and Th2 cells. T-helper cells that have not yet been presented to the antigen are denoted by 
ThO. In the model, two distinct routes (equivalent to two populations of antigen presenting 
cells) govern the maturing of a ThO cell: Antigen presented by the first (second) route elicit a 
ThI (Th2) response. In order to take into account the competition between mature Th cells, 
the induction ThO — > ThI (Th2) is forbidden when the neighborhood of the ThO cell has an 
overall majority of Th2 (ThI) cells. At last, a mature cell dies (it is substituted by a ThO 
cell) if such a cell is not restimulated by the appropriate antigen within the time interval Nt- 
The model can display a spontaneous symmetry breaking as one varies the antigen density 
or the cutoff Nt, in agreement with experimental results A probabilistic version of the 
model was proposed by Tome and Drugowich de Felfcio j(j by allowing the death of cells ThI 
and Th2 at each time step with a probability r. Thus, in this modified version of the model, 
the lifetime Nt was substituted by a mean lifetime related to the probability r. Also, the 
development of T H cells into either T H 1 or T H 2 cells occurs with a probability that depends 
on the neghborhood of the ThO cell and on a parameter p, related to the antigen density. Such 
probabilistic version presents the intrinsic spontaneous symmetry breaking found in the Brass 
et al. automaton, besides being amenable to an analytical approach [6|. The critical behavior 
of the probabilistic model proposed in Ref. p was studied in a subsequent work by Ortega et 
al. 0. By using a finite-size scaling analysis, from Monte Carlo simulations of square lattices, 
Ortega et al. determined the ratios of exponents ^ and ^, as well as the critical point of the 



2 



ed 

0, 



model. Their results suggest that the probabilistic automaton, although not satisfying detai 
balance, belongs to the same universality class of the two-dimensional kinetic Ising model 
thus supporting the up-down conjecture, introduced by Grinstein et al. |8[- In a subsequent 
work, Tome and Drugowich P] obtained the dynamic critical exponent z of the probabilistic 
automaton in two dimensions, by performing short-time Monte Carlo simulations by studying 
the collapse of the fourth order Binder's cumulant for different lattice sizes. 

In the present work we have revisited the probabilistic automaton of Ref. jf|. We have 
re-obtained the dynamic critical exponent z of the automaton by using a recent technique that 
mixes different initial conditions 10]. Such technique is based upon the short-time critical 



dynamics introduced by Janssen et al. ll| . who showed that even far from equilibrium the 



short-time relaxation of the order parameter follows a universal scale form 

M(t) = m t e } (1) 

where M(t) is the order parameter at instant t (measured in Monte Carlo steps per Spin - 
MCS), rriQ = M(0) is a small initial (t = 0) value of the order parameter, and 9 is the dynamic 
critical exponent, related to the increasing of the order parameter after the quenching of the 
system (When dealing with Monte Carlo simulations, the time t is discrete and measured in 
Monte Carlo Steps per Spin - MCS). Eq. (0) demands working with sharply prepared initial 
states with a precise value of mo- After obtaining the critical exponent 6 for a number of mo 
values, the final value for 6 is obtained from the limit m — > 0. 

Starting from an ordered state (m = 1), the order parameter M(t) decays in time, at the 
critical temperature, according to the power law [l^ 

(M(t)) mo=1 ~ (2) 

where ((•••)) i s the average of the quantity (•••) over different samples with initial order 
parameter value mo = M(0), (3 and v are the usual static critical exponents, related to the order 
parameter and to the correlation length, respectively, and z is the dynamic critical exponent, 
defined as r ~ £*, where r and £ are time and spatial correlation lengths, respectively. 

Starting from a disordered configuration with mo = 0, the second moment of the order 
parameter increases after the power law 

(M 2 (t)) mo=0 - t", (3) 

where the exponent to is given by 



and d is the dimension of the system. 

By combining Eqs. (J2J) and (J3J), da Silva et al. |l0| obtained the ratio 

which corresponds to a function with mixed initial conditions. At this point, it is important to 
stress here that, from Eq. (jSj) above, two independent runs are necessary in order to calculate 
the ratio F2. In one of them tuq = (numerator), while in the other one mo = 1 (denominator). 
The ratio F 2 has proven to be useful in determining the exponent z, according to recent studies 
of the 2D Ising model q = 3 and q = 4 states Potts models ^(J, Ising model with 
next-nearest-neighbor interactions jljj, Baxter- Wu model Q], at the tricritical point of the 
2D Blume-Capel model 2], at the Lifshitz point of the 3D ANNNI model Q| and in other 
studies concerning models with one absorbent state 17 1. 

In addition, in this work we have also calculated the global persistence exponent 9 g [18], 
related to the probability that the order parameter of the model, M(t), does not change its 
sign up to time t, after a quench of the system to the critical temperature. 

The layout of this paper is as follows: In section |H] we explain the model. In section ITTT1 we 
define the order parameter M(t) and we describe the methodology used in order to obtain the 
dynamic critical exponents 9 g and z. In section HVI the results obtained for the exponents z and 
9 g are shown. Finally, in section |S we present the main concluding remarks. 



II. THE MODEL 



In this work we have studied a two-dimensional probabilistic cellular automaton in which 
dynamics is governed by local stochastic rules. At each site i of the square lattice we have 
attached a variable <jj assuming the value 0, +1 or —1, depending on whether the site is occupied 
by a ThO, a ThI or a Th2 cell, respectively. Considering iV the total number of sites in the 
lattice, we have defined the set a = (cri, cr 2 , • ■ ■ , a^) to represent the microscopic state of the 
system. 

The probability of state a at time t, Pt{cr), evolves in time according to the equation 

P t+1 (a')=J2w(a'\a)P t (a), (6) 

a 

where the transition probability V^a'lcr) from state a to state a f obeys the condition 

£V(o» = l. (7) 

a' 
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On the other hand, for a system that evolves at discrete time steps, in which all the sites are 
updated at once, as is the case for the cellular automaton in this work, the transition probability 
VT(cr'|cr) is written in the form 

N 

W(a'\a) = l[u i (a! i \a), (8) 

i=i 

where uo^a'^a) is the conditional probability that site i be in the state a[ at time t + 1, given 
that the state of the system is a at instant t. Such conditional probability satisfies the condition 

= 1, (9) 

a' 

what implies immediately that Eq. (JZj) is fulfilled. The cellular automaton investigated in this 
work belongs to the class of totalistic cellular automata [19]. Thus, the transition probability 
LOi{a'^\a) depends on <Tj and on the sum £ = J2s &i+s, where the sum runs over the neighborhood 
of site i. More specifically, we are considering a particular kind of totalistic automaton for which 
uii^a'^a) depends only on the sign of the sum £. By defining 



Si = sign(£) = < 



1 if £ > 0, 

if £ = 0, (10) 
-1 if £ < 0, 



we may use the notation (^(cr^cxi, Si) in order to explicit the transition probability dependence 
on o"j and Sj. 

The transition probabilities (dynamical rules) are given by 

Wi(+1K Si ) = p5{ai, 0){5{ Si , +1) + -S(s h 0)} + (1 - r)5(a u +1), (11) 

Ui{-l\a h Si ) = P 5(a h 0){5( Si , -1) + 0)} + (1 - r)5(a h -1), (12) 

Wi(O|0i, Si) = (1 - p)5(a h 0) + r{5(a h +1) + 5(a h -1)}, (13) 

where, as discussed in section HI r is the death probability of T H 1 and T H 2 cells, and p is a 
parameter related to the antigen density 0. 

The dynamical rules (jll|) . (|12j) and (|13l ) above have "up-down", i.e., 

LUi(a'i\(Ji, Si) = uj^-a-l - <7i, -Si). (14) 



Therefore, following Grinstein et al. [8], we expect that the probabilistic cellular automaton 
investigated in this work be at the same universality class of kinetic Ising models. 
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III. THE ORDER PARAMETER AND THE DYNAMIC EXPONENT g 



A. Order parameter 

In the short-time Monte Carlo simulations performed in this work the order parameter M(t), 
as well as its higher moments M^ k '(t), were obtained from averages over a certain number of 
samples, N s . By defining the kth momentum of the order parameter in sample number j at 
instant t as 



(l 2 \ 



(15) 



the order parameter is written in the form 

N 3 /I? \ k 

M«(t) = <M*(t)) = JTTIk £ ( X>«W • ( 16 ) 

s j=l \i=l / 



As defined in Eqs. (J15|) and (jTSj) above, for k = 1 the order parameter is exactly the mean 
magnetization, M { - l \t) = (M(t)). 

Eq. (|T6|) defines the order parameter and its higher moments for a set of N s samples. Thus, if 
Nf, sets of samples are considered at instant t, there are Nb measurements of the magnetization 
(and its higher moments), M.l k \t), where 1 < I < N^. By considering such sets of samples, the 
final value of the magnetization and its higher moments are obtained from the average over the 
Nb sets of samples, i.e., 

N b 

MW(t) = (l/iV b )^M/ fc) (t), (17) 



i=i 



where Mf \t) is given by Eq. (|TH|) and the corresponding standard deviation is given by 

a M, W (t) = - V ( M, (fc) (t) - MW(t) ) 

L ' W J ^N b (N b - l) 2 V ' l V 

B. Dynamic critical exponent 9 g 

In this section we define the dynamic critical exponent 6 g and we describe two methods used 
in this work to obtain estimates for 9 g . 

In the first method , we have performed short-time Monte Carlo simulations in order to cal- 
culate the global persistence probability P(t), i.e., the probability that the magnetization does 
not change its sign up to time t. On the other hand, the probability P(t) is numerically equal 
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to the complement of the accumulated distribution p(t), according to which the magnetization 
changes its sign for the first time exactly at instant t, i.e., 

p® = i-EKO = i-E^, (19) 

t'=l f=l s 

where n(f) is the number of samples for which the magnetization changes its sign for the first 
time at instant t' and N s is the total number of samples. The exponent 9 g may be obtained 
directly from the power law scale relation [18] 

P(t)~r es , (20) 

from which we obtain lnP(t) = c — 9 g hit, where c is constant and each run requires a sharply 
prepared initial state, with a precise small value of mo, as discussed in Eq. (jTj). 

In a second way to obtain the exponent 9 g we have used the fact that the initial magnetization 
dependence of P(t) can be cast in the following finite-size scaling relation 

P(t) = r e * f{t/U) = L- 6 ^ f{t/L z ) , (21) 

which renders a different method to obtain the exponent 9 g from lattice sizes Li and L 2 18]. 
For this end we define W(t, L) = L egZ P(t), which turns out to be a function of t/L z . Therefore, 
if we fix the dynamic exponent z, the exponent 8 g can be obtained by collapsing the time series 
W(t 2 ,L 2 ) = f(t 2 /L 2 ) onto W{tx,Lx) = f(ti/Lf) as follows. Under re-scaling, with b = L 2 /L\, 
(L 2 > Li), we obtain 

W(t 2 ,L 2 ) = W(bH 1 ,bL 1 ) , (22) 
and the best estimate for 8 g corresponds to the minimization of 

x '(fl,) = Ef W(t ' L) -"T M) V (23) 
, \\W(t,L)\ + W(b'-t,bL) J 

by interpolating W to the time values b z t. In order to obtain the exponent 9 g using the collapse 
method described above, it is not necessary to fix a precise value of the initial magnetization m 
in the short-time simulations, once the scaling relation in Eq. (J2~T|) does not take into account 
the initial conditions of the system. So, we have used initial states in which (mo) ~ 0. On the 
other hand, the collapse method demands the dynamic exponent z to be known beforehand. 
Therefore, in this work we have used the scaling relation of Eq. (J5J) in order to obtain estimates 
for the exponent z. Although both methods described by Eqs. (|2()jl and (12 lj) were proposed in 
order to calculate the exponent 9 g of the two-dimensional Ising model (l8J, such methods were 
used recently for estimates of 9 g along the critical line and at the tricritical point of the 2D 
Blume-Capel model j^j. 
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IV. RESULTS FROM SHORT-TIME MONTE CARLO SIMULATIONS 



In this section we present details about the short-time Monte Carlo simulations performed 
for the cellular automaton considered in this work, as well as the results obtained for both 
dynamic critical exponents z and 9 g . 

A. Critical parameters 

Initially we refine the critical parameter r = 0.196 obtained in Ref. jfj) for p = 0.3. From 
Eq. (j2J we expected a straight line for the log-log plot of M(t) against t at the critical point 
(p = 0.3; r = 0.196). However, from log-log plots of Eq. (J2J) obtained for different values of r, 
we observed that more accurate straight lines are obtained for r ^ 0.196, as depicted in Fig. ^ 
for r = 0.190, 0.192, 0.194, 0.196 and 0.198. In the short-time simulations performed in order to 
obtain the curves shown in Fig d we have used square lattices (d = 2) with linear dimensions 
L = 160, N s = 10000 samples, iV 6 = 5 sets of samples and 1000 MCS. For each value of r used 




FIG. 1: Log-log plots of the order parameter (magnetization) M(t) versus t, constructed for p = 0.3 
and different values of r. The best linear fitting yields the best estimate for the critical value of the 
parameter r. From the results summarized in Tabled (see text), we have located the critical point 
(p = 0.3; r = 0.194). 

in Fig. ^we calculated the goodness of fit Q for the same time interval [t\, t 2 ] and obtained the 
values shown in Tabled from which we obtain the critical value r = 0.194. 

In order to confirm the critical value obtained in Table HJ we have also calculated the critical 
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TABLE I: Values of Q for different values of r obtained from linear fitting of log-log plots of the 
magnetization M(t) versus t. 



Time interval 


r 


Q 


[50, 300] 


0.190 


2.31 x 10~ 7 


[50, 300] 


0.192 


0.166 


[50, 600] 


0.194 


0.468 


[50, 300] 


0.196 


1.04 x 10~ 28 


[50, 300] 


0.198 


1.09 x 10~ 76 



value of r for p = 0.3 by using the effective exponent, which is given by 21] 

Such exponent takes into account finite-time corrections (finite-time scaling) and, in the limit 
t -> oo, Eq. fl23| behaves as H 

\{t) = ci + c 2 /t, (25) 

where c\ and c 2 are numerical constants and At is a fixed time step. 

Log-log plots of Eq. ()24|) are shown in Fig. for At = 10, p = 0.3 and r = 
0.190,0.192,0.194,0.196 and 0.198. From Fig. Hit is clear that the asymptotic behavior of 
X(t) given by Eq. (|25|) is verified for r = 0.194, thus confirming the result previously obtained 
from log- log plots of Eq. ©. 

Before presenting the estimates obtained for the critical exponents z and 9 g in the next 
subsections, it is important to emphasize here that the Monte Carlo simulations were performed 
only at the critical point (p = 0.3; r = 0.194). 



B. Critical dynamic exponent z 

In order to obtain the dynamic critical exponent z we have perfomed short-time Monte Carlo 
simulations for square lattices (d = 2) with linear dimensions L = 160, N s = 20000 samples, 
Nb = 5 sets of samples and 200 MCS. From the definition of the ratio F 2 given by Eq. (JSJ), 
we have constructed the log- log plot of F 2 {t) versus t shown in Fig. El obtained from Monte 
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FIG. 2: Effective exponent X(t) obtained for p = 0.3 and different values of r. From this graph, we 
have obtained the critical point (p = 0.3; r = 0.194) (See text for details). 

Carlo simulations performed for one set of samples. The straight line depicted in Fig. |3] is in 
accordance with the linear behavior predicted in the scaling relation of Eq. (jo]). 




FIG. 3: Typical log-log plot of ^(i) versus t (straight line), obtained from short-time Monte Carlo 
simulations at the critical point (p = 0.3; r = 0.194). 

In Table ITD we summarize the estimates for z obtained from different time intervals [ti,^]; 
together with the corresponding values of Q. 

From Table |H1 we have obtained the best estimate for the dynamic critical exponent z = 
2.097(8), once the goodness of fit Q — 0.99 is maximum at the corresponding time interval. 
However, such value of z is somewhat smaller than z = 2.17(2), which was obtained from the 
fourth order Binder's cumulant jjj]. In order to check the value of z obtained in this work, in 
the next section we obtain the global persistence exponent from the collapse method, which 
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TABLE II: Dynamic critical exponent z and the corresponding value of Q obtained for different time 
intervals. 



Time interval 


z 


Q 


[10, 120] 


2.0853(7) 


0.2 


[10, 160] 


2.0856(5) 


0.49 


[10, 200] 


2.0869(4) 


0.54 


[30, 200] 


2.084(4) 


0.98 


[100, 200] 


2.097(8) 


0.99 



depends on the value of z, and directly from Eq. (|20j) . which does not depend on the value of 
z. As we shall see in the following, the results for 9 g obtained from these two approaches are 
in very good agreement with each other. 

C. Critical dynamic exponent g 

By using Eqs. (|21|). (J22)) and (|23|) given in section 1111 Bl we have performed short-time 
Monte Carlo simulations in order to apply the collapse method of section 1111 Bl and obtain 
the dynamic critical exponent 9 g . According with the description of the method, the dynamic 
exponent z is to be known beforehand. Thus, we have used the value z = 2.097 obtained in 
the previous section. Monte Carlo runs were made up to 1000 MCS for N s = 40000 samples in 
square lattices with linear dimensions 20, 40 and 80. The collapse method was applied for pairs 
of linear dimensions (Li,L2) = (20,40) and (Li,L2) = (40,80), from which we have obtained 
6 g = 0.24(2) and 9 g = 0.25(2), respectively. In Fig. 0] we show the collapse of the curves 
obtained for L = 40 and L = 80 with 9 g = 0.25. 

Finally, we have obtained the dynamic exponent 9 g directly from the power law predicted 
in Eq. (J20|) . To this end, we have performed short-time Monte Carlo simulations to obtain the 
quantity P{t), from which we have constructed log-log curves of P(t) versus t, as depicted in 
Fig. E| From Eq. (j2"Uj) . the exponent 9 g was obtained directly from the slope of such curves. 
Monte Carlo simulations were performed for square lattices with linear dimension L = 80 and 
N s = 40000 samples, where each sample began from the initial state itlq = 0.005. Error bars 
were estimated from the Monte Carlo runs performed for Nj, = 5 sets of samples. In Table IIHI 
we present the results obtained for the exponent 9 g and the corresponding goodness of fit Q for 
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FIG. 4: Collapse obtained from the scaling relation Eq. (|21|) . corresponding to square lattices with 
linear dimensions {L\,Li = 40,80) and 6 g = 0.25. 




FIG. 5: Typical decay of the probability P(t) (straight line in a log-log plot), obtained from short-time 
Monte Carlo simulations of N s = 40000 samples with initial state mo = 0.005. 

three different time intervals. From Table ITTT1 ab ove . the best estimate for the global persistence 
exponent is 9 g = 0.247(4), corresponding to Q = 0.97. This result is in very good agreement 
with the estimates 9 g = 0.24(2) and 9 g = 0.25(2) obtained from the collapse method. 



V. CONCLUSIONS 



In this work we have revisited the probabilistic cellular automata proposed by Tome and 
Drugowich p] on the basis of a previous cellular automaton devised by Brass et al. Q]. From 
short-time Monte Carlo simulations, we have obtained the dynamic critical exponent z = 
2.097(8) by using a recent technique that mixes different initial conditions 10]. Such result is 
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TABLE III: Dynamic critical exponent g and the corresponding value of Q obtained for different time 
intervals. 



Time interval 


o 9 


Q 


[90, 300] 


0.247(4) 


0.97 


[90, 400] 


0.242(3) 


0.86 


[100, 500] 


0.238(6) 


0.95 



slightly smaller than the value z = 2.17(2) obtained from the fourth order Binder's cumulant 

0. 

We have also performed short-time Monte Carlo simulations in order to obtain estimates 
for the dynamic critical exponent 6 g , the global persistence exponent, by using two distinct 
approaches: The collapse method described by Eqs. (J21j) . (|22j) and (|23j) . and directly from the 
power law scaling given by Eq. (|2U|) . From the collapse method, by using square lattices with 
linear dimensions (Li,L 2 ) = (20,40) and (L 1 ,L 2 ) = (40,80), we have obtained 9 g = 0.24(2) 
and 9 g = 0.25(2), respectively. Directly from Eq. (J20j) we have obtained 6 g = 0.247(4), in very 
good agreement with the estimates yielded from the collapse method. 
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